Longitudinal Survey 3 Residual Pools

In this script we estimate residual pool characteristics. The residual pool is the pool that would form if streamflow were to decrease to zero. Its depth is characterised as the elevation difference between the pool’s deepest point and the pool outlet / riffle crest that controls the flow of water out of the pool. It serves as a flow-independent metric of pools in a stream and, with several assumptions, we can estimate the depths, areas, and volumes of residual pools based on our field data. See Lisle (1989) for more background on the concept of residual depths.

We tried to make field data collection as simple as possible by relying on water depths rather than elevations for our longitudinal profiles, so that people could collect data without survey equipment - reducing the cost and expertise required, and facilitating assessment in densely vegetated / woody coastal streams, highly incised channels etc. The use of simple depth data requires a couple of extra processing steps compared with using elevation data, but luckily some rather clever individuals have already figured it out.

The code we use here is adapted from the ‘Rapid Bed Profile’ approach from Stack (1989), with additional guidance from Kaufman et al. (1999), and we adopt Stack’s slope factor. The slope factor was developed in coastal streams but has been applied successfully in a variety of other systems. If you have the resources, the accuracy of the slope factor for your particular system can be tested by comparing water-depth-derived to elevation-derived residual pool data. However, even if this assumption is untested (as it will be in most cases), the resulting metrics remain valid for assessing changes within a system, as long as you don’t change the slope factor between years or reaches in your comparison.

Data import

First we will import the tidied dataframe generated by the Data Wrangling script, which was exported to your output folder. This file has longitudinal thalweg water depth data that includes both regular predetermined spacing and any ‘extra’ pool outlets and max pool depths encountered (‘XPO’, ‘XD’). We proposed that the ‘extra’ measurements were taken in the field, in order to capture the depth minima (potential pool outlets) and maxima (deepest pool sections) that might occur between regular intervals.

In this code we will require some additional constants that we ask you to define in the code below based on data you recorded in the field (or are able to estimate). This is done by creating a dataframe/table, for which you must enter the following:

• names of each channel
• the years data were collected
• the standard interval spacing (in metres) for each reach and each year
• the average reach slope (0.5 = 0.5%, not 50%) for each reach and each year
• the average wetted width during low flow (in metres) for each reach and each year.

Enter this information in the code below. Data is entered in vectors of equal length and following the same sequence. For example, in the existing code, Hatchery Channel is associated with spacing of 1 and 1, slope of 1 and 1.1, and wetted width of 6.8 and 7.0, for 2024 and 2025 respectively. Once you create the dataframe with your data, check through the table to ensure these constants are assigned to the correct reach & year.

# Create a data frame to hold constants for multiple sites and years. Enter your data below:
constants_df <- data.frame(
  Site = c("Hatchery Channel", "Hatchery Channel", "Brousseau Channel", "Brousseau Channel"),  # Edit / Add your site names
  year = c(2024, 2025, 2024, 2025),  # Corresponding years for each site
  spacing = c(1, 1, 1, 1),  # spacing values for each site-year combination
  Slope = c(0.5, 0.6, 0.5, 0.5),    # Slope values for each site-year combination, in %
  WetWidth = c(6.8, 7.0, 6.5, 7.1)  # Average wetted width in m, for each site-year combination
)

# View the constants data frame
print(constants_df)

##                Site year spacing Slope WetWidth
## 1  Hatchery Channel 2024       1   0.5      6.8
## 2  Hatchery Channel 2025       1   0.6      7.0
## 3 Brousseau Channel 2024       1   0.5      6.5
## 4 Brousseau Channel 2025       1   0.5      7.1

Creation of Residual Surfaces

The residual surface is the water surface at the point when streamflow approaches zero, and residual depth is the vertical distance between the residual surface and the stream bed. As such, when we estimate residual pools and residual depths, these are features that are independent of flow conditions at the time of the field survey. While they might be representative of actual pool conditions during times of very low flow, their value lies with the ability to make comparisons of pool habitat among years and/or reaches without quantifying and controlling for variable flow conditions.

Our code uses the sequence of thalweg depth measurements to identify potential pool outlets in the sampled reach, and estimates residual depths at each measurement location based on a projected residual surface that extends upstream from each pool outlet. The residual surface is projected using the expected pool water surface slope (which is estimated using Stack’s empirical equation, incorporating the average reach slope that you entered as a constant).

It is important to note that when we identify pool outlets or riffle crests (either in the field or using this code), we usually base this determination on the bed elevation and/or water depth at the current location compared with the areas immediately adjacent. However, a potential pool outlet will not ‘behave’ as a pool outlet under all flow conditions: it could be high and dry under very low flow conditions, or it could itself be within a pool due to the backwatering effects of a higher-elevation pool outlet downstream. From here onward, we will use the term ‘local depth minima’ (‘LDM’ in code) to refer to locations that have the potential to be pool outlets - that is, they are less deep than areas immediately downstream and upstream, accounting for the slope of the stream. We will reserve the term ‘pool outlet’ for cases when local depth minima would be expected to actually backwater a pool under flow conditions approaching zero.

At this point our code will generate a new dataframe (residuals.df) with the local depth minima identified and the residual surface projected. If you are interested in how this is done, check out the R code. For everyone else, here are some fish:

🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟 🐟🐟🐟🐟🐟

Below we generate a quick visual of the depths and residual surfaces that were projected by the code. Take a look to ensure it matches your understanding of residual pools and look for potential errors (the figure is interactive, so zoom in if you have several plots). Try to identify some local depth minima that are within a pool (if present in your data), and compare these to the local depth minima that act as pool outlets - the shallower features at which residual surfaces begin. Remember that the residual surface slopes down in these figures as you progress upstream (right in the plots), which is to compensate for the fact that the true thalweg surface slopes down as you progress downstream.

‘Meaningful’ Pools

Consider the vertical exaggeration of the above figure - some ‘pools’ may be very shallow indeed, and might not necessarily correspond with what you would actually consider a pool in a stream.

It is probably most informative to consider metrics relating only to those pools that are expected to be hydraulically and/or biologically meaningful. The code now detects pools and assigns each a pool ID, so that we can characterise them individually and select subsets as appropriate.

The above figure should show a series of distinct pools and is exported to your output folder. Don’t worry if the shading cuts off before the pool’s upstream extent, this is just an artefact associated with my limited figure-plotting abilities. The important thing to note is that the calculations below will integrate the full area beneath the residual surface for each pool.

Pool Metrics

The code now calculates the following metrics for each residual pool:

• Maximum depth;

• Pool length (pool outlet to where the residual surface meets the thalweg at the upstream end of the pool); and

• Sagittal area.

Sagittal refers to a plane that divides a body into left and right, so imagine a two-dimensional vertical screen extending downstream all along the thalweg, separating river left from river right. Within each residual pool, the sagittal area corresponds to the surface area of that screen. Too abstract? Take a look at the above figure again: Each of the shaded areas corresponds to the sagittal area for that particular residual pool - sagittal area integrates each depth measurement over the length of the pool. Although this may not be as intuitive as pool volume, it requires fewer assumptions than will be necessary to generate estimates of pool volume (i.e., regarding cross-sectional channel shape of pools at zero flow, wetted widths). Because fewer assumptions are required, and because it easily corresponds to 2D visualisations like the one above, we recommend using sagittal area rather than volume for a metric tracking changes in pool habitat (though we will also calculate volume later).

Subsetting Pools

So far we have identified all of the residual pools, but depending on your dataset some of the pools visualised above might look more like puddles. So to focus our summary statistics on pools that better align with what we consider biologically- and/or geomorphically-meaningful, we can use the metrics we just calculated to focus on a specific subset of deeper pools.

In terms of subsetting pools, some people use a specific depth cut-off for a species of interest (e.g., Mossop and Bradford 2006 used 0.1 m). Alternatively the pools that we focus on can be based on physical characteristics of the sampled reach. We will do the latter, keeping this code more widely applicable, by following Stack (1989) who determined meaningful pools as those that had a maximum depth greater than or equal to the standard deviation of depth for the reach.

In the following code chunk it reads ‘custom_depth <- NULL’. This is our default subsetting (selects only pools with a maximum residual depth greater than the standard deviation of the residual depths for the entire sampled reach). If you want to subset the pools in another way, you can modify the code here. For example, if it is biologically meaningful for your system to consider only pools that are 0.1 m deep or greater during approaching-zero-flow periods, you can specify ‘custom_depth <- 0.1’. Or, if you want to consider all pools no matter how shallow, then assign custom_depth a value of ‘0’. Just be sure to compare like-for-like when you consider different years/sites that you may wish to compare.

# Add a custom depth threshold in metres, or leave as 'NULL' to use standard deviation. If you are interested in all of the residual depths (i.e., not subsetting), just enter '0' for the custom_depth
custom_depth <- NULL

The table above has been exported to your output folder. It displays the subset selection of residual pools, the sagittal area of each pool, the maximum residual depth (deepest point per pool at zero flow), and the pool length of each pool. When considered at the reach-scale, summary statistics of the maximum residual depth and pool length have been considered useful metrics to track pool quality and quantity (e.g., Lisle 1986, Mossop and Bradford 2006, Clark et al 2019).

The reach mean of the pool maximum residual depths has been used as a proxy for pool quality, and positive relationships with Chinook salmon density have been demonstrated (e.g., Mossop and Bradford 2006). Relationships between habitat metrics and biota can be highly context-dependent, and independent verification is recommended where possible. However, as the resources required to test these assumptions are often prohibitive, we recommend taking advantage of relationships that have already been quantitatively demonstrated in the scientific literature. Although it is possible to generate quantitative estimates of biotic densities if the relationship between habitat metric and biotic density is monotonic, we feel that this is too speculative without data for your individual system. However, for the purpose of restoration monitoring, we feel it is incredibly powerful to be able to demonstrate (e.g.) that we have caused an x-sized increase in a habitat metric that is known to be positively correlated with target species density in a similar system.

The proportional length of pools in the reach also has demonstrated positive relationships with coho, Chinook densities (Clark et al 2019). Bear in mind that there will always be nuance involved in the application and interpretation of habitat metrics. You must understand and predict, in an a priori fashion, what your restoration objectives are with respect to habitat: Does your reach need more pool of a certain depth? Is your reach lacking much-needed riffles? Presumably there are cases where the reach has too much pool habitat. These types of questions are best addressed with well-matched reference sites, but local knowledge and the scientific literature will also be of great help.

We will now generate and export to your output folder some summary statistics, including the aforementioned mean max residual depth and length in residual pool, for your data:

These above statistics should be useful for tracking changes over time, or comparing reference / control reaches. Because they are all based on residual depths, they are flow-independent.

Total sagittal area and proportion of the reach length in residual pool should be compared between years and/or sites, and you can interpret whether the changes are in the predicted direction and of a meaningful magnitude. Since these two metrics have one value per reach, the metrics themselves cannot be evaluated statistically (unless you have a large number of reaches/years), however a directional change over time that aligns with your predictions is obviously quite a good sign.

Inferential Tests

We can apply inferential statistical tests to the sagittal area, maximum pool depth, and pool length data - that is, the data which have values at the level of each pool (i.e., a population of values is available per reach/year). We will use inferential statistical tests to determine whether the differences we observed among reaches and years are likely to be due to chance. You will be familiar with p values of 0.05 (observed result expected just due to chance 1 in 20 times) and we will use p=0.05 as our cut-off here. However, depending on your adaptive management plans, you may be willing to accept different levels of certainty, and p values of 0.1 are not unheard of in restoration projects. This all depends on the risk tolerance of you and your team, and should be established before conducting analyses.

For these inferential tests we will use ANOVAs because of their familiarity for ecologists and expected suitability for the data. For the metrics relating to residual pools that we have generated, we will test for differences between sites, years, and an interaction of sites*years (i.e., did the different locations change differently over time). The code uses the subset of pools that were defined earlier as ‘meaningful’.

Choice of Metric: At this point, you should choose which of ‘sagittal area’, ‘maximum pool depth’, and ‘pool length’ best aligns with your needs (e.g., relevance to restoration objective, comparisons with literature, ease of communication). You will base your interpretation of the ANOVA on this chosen metric, even though we will run ANOVA for all three.

The reason for selecting one metric for interpreting your hypothesis is that the metrics are not independent of one another: Sagittal area incorporates depth and length components, and all three metrics are connected in the stream; if you used a subset of the residual pools, all metrics inherently include some information about depth. If you were to continue without deciding upon a metric, it would be too tempting to come up with post-hoc explanations (e.g., if sagittal area increased but not depth or length). You should certainly consider potential implications of your results (carefully), but in terms of evaluating the success of your project, it is important that you decide in advance which metric change = success.

Although it may be possible to run a MANOVA test (which can accommodate correlated dependent variables) using all three metrics together, MANOVA tests typically require larger sample sizes than we anticipate for many reach surveys. Therefore, we stick with multiple ANOVAs with a large dose of caution.

Assumption Robustness

Whenever applying statistical tests, we should ensure that our data is suitable, in that it meets the assumptions of the intended inferential test. For ANOVAs the assumptions are:

• Data in different groups should be independent:

  • One of your sites must not cause changes in another (though both sites experiencing shared ‘external’ causal factors, such as being in the same watershed, is acceptable). This is hopefully accounted for by good selection of impact/control/reference sites.

  • You may note that no site is actually independent from that same site in the past. For longer time-series we would need to adopt a different inferential approach to avoid autocorrelation, but here we are dealing with simple before vs after comparisons, and it is widely accepted that ANOVAs are suitable for before/after comparisons if the comparison is balanced (e.g., 2 years of before data compared with 2 years of after data).

• Residuals are normally distributed: In this case we mean residuals more broadly - the difference between each datapoint and the mean. We will examine this assumption using Q-Q plots and the Shapiro-Wilk normality test.

• Homoscedasticity: Variance should be broadly equal across all groups. We will test this with Levene’s test.

If any of these assumptions are violated, we can transform the data and test whether the transformed data (log-, arcsine- etc.) satisfies the assumptions. If this is not an option, we should seek a more suitable statistical test (e.g. Welch’s ANOVA or mixed-effects models).

The tests of assumptions are carried out below, after we specify the before and after data that you wish to compare.

Before / After Years

Enter your before and after years in the below code, you may have more than one year in each category. Remember that these should balance (equal number of years before as after), but they need not be consecutive years. If you do not have balanced data (e.g., only one year was possible before project implementation) be very cautious in running and interpreting an ANOVA. You could select only a subset of your data to keep the comparison balanced, but if you do this for more than one set of dates you must account for the multiple comparisons inflating the likelihood of a Type I error (false positive).

# Define "Before" and "After" years
before_years <- c(2024) #add more years inside brackets if needed, e.g. 'c(2024, 2025, 2026)'
after_years <- c(2025) #add more years inside brackets if needed, e.g. 'c(2027, 2029, 2036)'

The next section of code generates tests and figures that you should review to ensure your data are suitable for an ANOVA. In brief, look for the following:

• Boxplots: consider whether outliers (dots) are real or errors
• Q-Q plots: a reasonably straight line of points means alignment with normal distribution
• Normality test (Shapiro-Wilk): p<0.05 indicates non-normality
• Levene’s test: p<0.05 indicates that variance among groups differs
• Skewness values < 0 are left skewed, > 0 right skewed. Kurtosis values < 3 light-tailed, > 3 heavy-tailed.

Normality Test Results for Pool Summary Data

Site	Time Category	Max Residual Depth P-Value	Sagittal Area P-Value	Pool Length P-Value	Message
Brousseau Channel	After	0.700	0.346	0.507	Normality test completed
Hatchery Channel	After	0.430	0.325	0.977	Normality test completed
Brousseau Channel	Before	0.471	0.129	0.086	Normality test completed
Hatchery Channel	Before	0.172	0.081	0.051	Normality test completed

## 
## --- Levene's Test Results for sagittal_area ---
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  3  6.2061 0.002838 **
##       24                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## --- Levene's Test Results for max_residual_depth ---
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.9027 0.05559 .
##       24                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## --- Levene's Test Results for pool_length ---
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.7929 0.1754
##       24

##                              Variable  Skewness Kurtosis
## sagittal_area           sagittal_area 1.7710682 4.904845
## max_residual_depth max_residual_depth 0.8824034 2.069536
## pool_length               pool_length 1.3271828 4.412251

If, based on the above evaluation, your data appear to be unsuitable for the intended inferential test you may wish to use an alternative non-parametric test or transform the data. We provide some assistance for transformations here, but this part may require some more work on your part since we cannot predict what your data will need.

As an example, our data appeared normally distributed (p>0.05 Shapiro-Wilk) but Levene’s test revealed that sagittal area variance was unequal across site-years (p<0.05). We therefore log-transformed the sagittal area data and re-ran the above diagnostics, with Levene’s test returned as non-significant and therefore acceptable to be used in ANOVAs. Note that only those metrics that do not meet the assumptions need to be transformed.

The code we used for the transformation is below. If you need to transform the data, follow these steps:

Transformations

1. Add your transformation(s) to the code chunk below. You must add a formula to calculate an appropriate transformation for each applicable variable within transformations <- list(…)). As an example, if you want to try a log transformation, replace ‘identity,’ in the below code with ‘function(x) ifelse(x > 0, log(x), NA),’

2. Re-run the above diagnostics code after updating the applicable variable name(s) for each diagnostic test: where it says ‘variables <- c(“sagittal_area”, “max_residual_depth”, “pool_length”)’ you will update each applicable transformed variable name(s) to “transformed_variable_name”

3. Once you have a transformed variable that meets the assumptions of the ANOVA, update the ANOVA inferential test code (in the subsequent code chunk) to use the “transformed_variable_name”

# Function to apply transformations and append to the dataframe
apply_transformations <- function(data, transformations) {
  for (var in names(transformations)) {
    transform_fn <- transformations[[var]]
    transformed_var <- paste0("transformed_", var)
    # Apply the transformation and append as a new column
    data[[transformed_var]] <- transform_fn(data[[var]])
  }
  return(data)
}

# Define transformations as a named list. You can add appropriate transformations here, then once you confirm the transformation is suitable, amend the anova code to call the correct variable

# Use 'identity' if no transformation is needed
transformations <- list(
  sagittal_area = identity, # example code fo log-transformation: replace 'identity,' with 'function(x) ifelse(x > 0, log(x), NA),' 
  max_residual_depth = identity,  # No transformation
  pool_length = identity # no transformation 
)

# Apply transformations and update the dataframe
pool_summary_filtered <- apply_transformations(pool_summary_filtered, transformations)

ANOVA

Once you have satisfied yourself that your data (or transformed data) are suitable, we will run the inferential tests. We use ANOVA to ask whether the following metrics differ among sites, years, and site*year interactions:

• Sagittal area
• Maximum residual depth
• Pool length

If you did use transformed data, ensure you reference that variable in the below code (e.g., change “sagittal_area” to “transformed_sagittal_area”)

# Function to run ANOVA for a given variable
run_anova <- function(data, variable) {
  formula <- as.formula(paste(variable, "~ Site * Time_Category"))
  anova_result <- aov(formula, data = data)
  summary(anova_result)
}

# Ensure transformed data is used in the ANOVA if necessary (change the variable name in quotation marks in these next rows)
anova_max_residual_depth <- run_anova(pool_summary_filtered, "max_residual_depth")
anova_pool_length <- run_anova(pool_summary_filtered, "pool_length")
anova_sagittal_area <- run_anova(pool_summary_filtered, "sagittal_area")

## [1] "_______________ANOVA for Max Residual Depth______________"

##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                1  0.072   0.072   2.907   0.1011    
## Time_Category       1  3.712   3.712 150.194 8.07e-12 ***
## Site:Time_Category  1  0.115   0.115   4.642   0.0415 *  
## Residuals          24  0.593   0.025                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "_______________ANOVA for Residual Pool Length______________"

##                    Df Sum Sq Mean Sq F value  Pr(>F)   
## Site                1      6     6.0   0.038 0.84625   
## Time_Category       1   2022  2021.5  12.868 0.00148 **
## Site:Time_Category  1    163   163.4   1.040 0.31792   
## Residuals          24   3770   157.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "_______________ANOVA for Residual Pool Sagittal Area______________"

##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                1   31.6    31.6   0.675    0.419    
## Time_Category       1 1388.2  1388.2  29.659 1.35e-05 ***
## Site:Time_Category  1    6.3     6.3   0.134    0.718    
## Residuals          24 1123.3    46.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Take some time to ensure you understand what the above tells you, particularly with regards to your chosen metric of success. In the example ANOVA above, there was one site:time interaction that was significant at p=0.05 (maximum residual depth), and the time category was significant for all metrics. This suggests that a change in metrics over time was not likely due to chance, but only for maximum residual depth did our restoration site change over time differently than the control site. However, we chose sagittal area as our a priori success metric, so we cannot interpret this as a successful restoration action at this point in time.

Although we could not discard our alternative hypothesis (that changes at our restoration site were not due to our actions), we may still learn from the above to inform future plans. Consider all the information we have generated so far for you data, and take some time to think about it. Step outside to do this if you are able to. Consider what it might mean if pools at our restoration site appear deeper on average, but sagittal areas did not change? Consider whether there are pools that did or did not scrape through the subsetting cut-off (i.e., very small pools of marginal value), how might they impact your mean values? Was there a big enough change in these marginal cases to create statistical significance? Is any of this biologically or geomorphically meaningful based on your hypotheses?

This is all endlessly interesting/maddening (depending on your perspective) and we highly recommend finding another person (or several) to talk it through with. Before moving on to further analyses, ensure you have a reasonable understanding of the information so far. Take a step back to think about the big picture.

Volume Estimation

If you have no hypotheses or predictions specifically regarding pool volumes, you may wish to skip this section and move onto Fine Sediment in Pools.

If you took wetted width measurements at transects along the reach, you can use the average to estimate the volume of residual pools. This approach is a more simplified approach than that of Robison (1998) who linked wetted width measurements to nearby longitudinal intervals. If this is possible for the data you collected, we recommend reviewing Robison (1998) to obtain more ‘localised’ estimates of wetted width. However, bear in mind that there may regardless be a bias in residual wetted widths if (e.g.) pools tend to be wider than riffles. For our simplified approach we assume that the wetted width is constant for the reach, and we also assume that the stream cross sections are perfect triangles. We then use the ratio of residual depths to field-recorded depths to estimate the residual wetted widths of pools when flow approaches zero.

Considering changes in (coarsely-estimated) pool volume can be of interest due to the potential for natural processes to alter both the longitudinal and cross-sectional heterogeneity of a stream. For example, perhaps we are anticipating a narrower wetted channel with deeper pools, and want to explore what this means for pool volume. However, as suggested by the coarse assumptions we made for this calculation, interpret results with caution.

The above table lists the volume estimates for each of the ‘meaningful’ subset of residual pools, and this table is also exported to your output folder. Next we will display and export a table with some key summary statistics for the reach: Total volume of residual pools, and the mean and standard deviation residual volume of pools.

Inferential Tests - Volume

Following the same approach as for the one- and two-dimensional pool metrics, we will use ANOVAs to test the pool volume differences between sites and years.

Assumption Robustness

As before, the first step is to check whether an ANOVA is suitable for the variables.

Normality Test Results for Pool Volume Data

Site	Time Category	Pool Volume P-Value	Message
Brousseau Channel	After	0.500	Normality test completed
Hatchery Channel	After	0.461	Normality test completed
Brousseau Channel	Before	0.451	Normality test completed
Hatchery Channel	Before	0.014	Normality test completed

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  3  6.9518 0.001577 **
##       24                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Skewness: 1.814195

## Kurtosis: 5.117599

We needed to transform our volume data, which was not normally distributed. Following the same transformation steps as before, We replaced ‘pool_volume_m3’ with the log-transformed ‘transformed_pool_volume_m3’ (using the below code) and re-ran the diagnostics. Our Q-Q plot and tests looked much better, so we used the transformed variable in the ANOVA.

# Define transformations as a named list for pool volume
transformations <- list(
  pool_volume_m3 = function(x) ifelse(x > 0, log(x), NA)  # Log transformation with handling for non-positive values (volumes should all be positive!)
)

# Apply transformations and update the dataframe
pool_volumes_filtered <- apply_transformations(pool_volumes_filtered, transformations)

ANOVA - volume

# Ensure transformed data is used in the ANOVA function if necessary (change the variable name in quotation marks)
anova_pool_volumes <- run_anova(pool_volumes_filtered, "pool_volume_m3") # adjust to "transformed_pool_volumes_m3" as needed

##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                1      4       4   0.014    0.907    
## Time_Category       1   9095    9095  33.011 6.39e-06 ***
## Site:Time_Category  1    231     231   0.839    0.369    
## Residuals          24   6613     276                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

When interpreting potential changes in pool volume metrics, remember that we needed to make assumptions regarding the widths of residual pools, making these volume metrics more speculative than sagittal area metrics. While volumes can be useful for interpretation, if there are disagreements between the inferential test results among volume- and area/length-based metrics, we recommend defaulting to the latter (unless, of course, you have a solid a priori hypothesis and prediction that suggests volume should differ in the way you observed).

Fine Sediment in Pools

The longitudinal survey included presence/absence observations of fine sediment along the thalweg. Here we will generate a metric of the proportion of residual pool thalweg that features fine sediments.

Note that, while the previous metrics in this script use the ‘custom depth’ threshold to subset what are considered ‘meaningful’ pools, we recommend considering fine sediment relative to the full contingent of residual pools that the code has identified. Our rationale is that the fines metrics relate more closely to the mechanisms that can be responsible for making pools ‘meaningful’ (or not) in terms of fish habitat availability. If we were to look at the proportion of fines only in ‘meaningful’ (deep enough) pools, we might miss pools that are closer to being ‘not pools’ (i.e., recently filling-in or scouring of sediment). Of course, when pools are fully infilled with sediment they will not be detected as pools at all, but the ability to recognise when a pool is filling in (or scouring out) may be an important indicator of restoration success, or a trigger for adaptive management actions.

However, should it be of greater interest for your project, we have included in the code below the ability to select a subset of pools. Consistent with the previous custom depth threshold, you can add a specified depth in metres, or enter ‘NULL’ to select only pools with a maximum residual depth greater than the standard deviation of the residual depths for the entire sampled reach. Our default is 0 (no subsetting).

# Define custom threshold if required. We recommend 0 (all residual pools, regardless of maximum depth)
custom_depth <- 0  # Replace with your desired threshold in metres, or set to NULL to use the standard deviation of the reach (that year)

Brousseau Channel - Pool Fine Sediment

Site	year	pool_id	pool_length	fines_length	fines_proportion
Brousseau Channel	2024	1	21.5	21.0	0.977
Brousseau Channel	2024	2	0.5	0.5	1.000
Brousseau Channel	2024	3	21.5	17.5	0.814
Brousseau Channel	2024	4	0.5	0.5	1.000
Brousseau Channel	2024	5	0.5	0.0	0.000
Brousseau Channel	2024	6	8.5	5.0	0.588
Brousseau Channel	2024	7	4.5	4.5	1.000
Brousseau Channel	2024	8	8.5	7.0	0.824
Brousseau Channel	2024	9	2.5	1.0	0.400
Brousseau Channel	2024	10	22.5	20.0	0.889
Brousseau Channel	2024	11	1.5	0.5	0.333
Brousseau Channel	2024	12	2.5	1.0	0.400
Brousseau Channel	2024	13	2.5	0.0	0.000
Brousseau Channel	2024	14	12.5	8.0	0.640
Brousseau Channel	2024	15	3.5	3.5	1.000
Brousseau Channel	2024	16	6.5	6.5	1.000
Brousseau Channel	2025	1	61.0	50.5	0.828
Brousseau Channel	2025	2	18.0	13.0	0.722
Brousseau Channel	2025	3	33.0	23.0	0.697
Brousseau Channel	2025	4	30.0	21.5	0.717
Brousseau Channel	2024	Total	120.0	96.5	0.804
Brousseau Channel	2025	Total	142.0	108.0	0.761

Hatchery Channel - Pool Fine Sediment

Site	year	pool_id	pool_length	fines_length	fines_proportion
Hatchery Channel	2024	1	24.5	0.0	0.000
Hatchery Channel	2024	2	13.5	0.0	0.000
Hatchery Channel	2024	3	10.5	0.0	0.000
Hatchery Channel	2024	4	5.5	0.0	0.000
Hatchery Channel	2024	5	30.5	29.0	0.951
Hatchery Channel	2024	6	11.5	11.5	1.000
Hatchery Channel	2024	7	12.5	12.5	1.000
Hatchery Channel	2024	8	8.5	8.5	1.000
Hatchery Channel	2024	9	9.5	9.5	1.000
Hatchery Channel	2024	10	1.5	0.0	0.000
Hatchery Channel	2025	1	39.0	0.0	0.000
Hatchery Channel	2025	2	19.0	0.0	0.000
Hatchery Channel	2025	3	1.0	0.0	0.000
Hatchery Channel	2025	4	57.0	55.5	0.974
Hatchery Channel	2025	5	24.0	20.5	0.854
Hatchery Channel	2024	Total	128.0	71.0	0.555
Hatchery Channel	2025	Total	140.0	76.0	0.543

The above tables (also exported to output folder) include a final row for each site-year with reach total pool length, reach total pool length with fine sediment, and the reachwide proportion of pools with fines (a ratio from 0 to 1). Comparison of the reachwide proportion alone may be sufficient to demonstrate project effectiveness, perhaps complemented with the figure that we will create below. However, if you would like to statistically compare how fine sediment changed over time compared with a control, we can look at the distributions of underlying pool data (the fines_proportion data for every pool).

Inferential Tests - Fines

This time, we are dealing with data that are bounded by 0 and 1. And in many cases, we also anticipate a relatively high frequency of 0s (no fines in pools) and/or 1s (fines throughout). We are confident that non-transformed data would violate the assumptions of ANOVA, so we omit the initial diagnostic tests. Our attempts to arcsine-transform the data were not successful, as our transformed data remained non-normally distributed with significant heteroscedasticity. As such, and since we suspect our fines data distribution is not particularly unusual, we opted for Beta regression.

Slightly nerdy bit, feel free to ignore: Beta regressions are capable of dealing with non-normal heteroscedastic data. Although they cannot actually handle 0’s and 1’s, for our case they are preferable over two-part models (e.g. hurdle models): Since an observation of “no fine sediment” could very plausibly be recorded where 2 cm of fine sediment was actually present in a 20 m pool, we are comfortable adjusting the 0’s and 1’s to 0.001 and 0.999. We feel a hurdle model would excessively complicate interpretability (e.g. a binomial fines presence/absence mechanism plus a fines quantity mechanism). See Geissinger et al 2022 if interested.

We will generate a frequency plot of the fine proportion data, fit the Beta regression, and plot the frequency distribution of the residuals (which should be broadly symmetrical around zero) and their relationships to the fitted values (for which no clear relationship or heteroscedasticity should be evident). We will only export the results of the Beta regression to your output folder.

## Warning in betareg.fit(X, Y, Z, weights, offset, link, link.phi, type, control,
## : no valid starting value for precision parameter found, using 1 instead

## 
## Call:
## betareg(formula = fines_proportion_adjusted ~ Site * Time_Category, data = fines_proportion_filtered)
## 
## Quantile residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6659 -0.7903  0.0613  1.0070  1.4205 
## 
## Coefficients (mean model with logit link):
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                0.2570     0.6770   0.380    0.704
## SiteHatchery Channel                      -0.9495     0.9034  -1.051    0.293
## Time_CategoryBefore                        0.1700     0.7544   0.225    0.822
## SiteHatchery Channel:Time_CategoryBefore   0.4271     1.0473   0.408    0.683
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)  0.52358    0.09579   5.466 4.61e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 54.34 on 5 Df
## Pseudo R-squared: 0.1048
## Number of iterations: 20 (BFGS) + 2 (Fisher scoring)

Since you may not be familiar with Beta regression, here are some basics on interpretation for the test output:

• “Optimization failed to converge…” If you see this warning, it may be that there were 0 values across an entire Site. This should be fairly easy to interpret without inferential tests.

• Quantile residuals measure the difference between the actual data and predicted values. The median should be close to 0. Also check the residuals vs fitted plot to ensure there is no clear trend or unusual outliers

• Coefficients (mean model): these are the effects of the predictors (site, time) on the response variable (fine sediment proportion in pools). Look for significant p values in the last column, particularly for the Site_____:Time_Category_____ row, which is the interaction term (if significant, this means the proportion of fines in residual pools changed over time differently between sites).

• Phi coefficients: there is a p-value here as well, but this refers to whether the spread of the data around the mean was an important feature. If the phi estimate is very high, this means the data was concentrated around the mean (fines proportion doesn’t fluctuate much between sites/times). Consider phi of 0.1 to be spread out data, 1 to be moderate, and 10 to be concentrated around the mean.

• Pseudo R-squared: this is how much of the variation in fines in pools can be explained by the site and time variables. Often low in ecological contexts.

Proportional data can be less intuitive to interpret, and using an unfamiliar inferential test might not clear things up. Not to worry: After we perform a similar analysis with the aquatic vegetation data, we will create some plots that show the distribution of fines and aquatics across the reach. These plots are more intuitive and will help you to interpret any changes in greater context.

Aquatic Vegetation in Pools

In the field, you may have also estimated the proportion of the wetted width that was occupied by aquatic vegetation. During data wrangling we added three new data columns to facilitate analyses of these aquatic vegetation observations: subveg (submerged vegetation), floatveg (floating vegetation), and emergveg (emergent vegetation). We will use these variables to estimate the length and plan area of each residual pool that would be occupied by aquatic vegetation when flows approach zero. These can be useful metrics for tracking habitat quality and mechanisms of habitat change, however plan area estimates come with a caveat:

In estimating plan areas of vegetation, we assume that the proportion of residual pool width occupied by vegetation is equal to the proportion of the wetted width that was occupied in the field. In most cases this will be an overestimate because the thalweg will typically have less vegetation and, as flows recede, the thalweg makes up a greater proportion of the wetted width. Although vegetated residual pool areas can be useful for tracking changes over time, because they are based on proportions of wetted widths, you should not consider them flow-independent. When interpreting changes over time (or from different reaches), remember that if flows were notably different, then the amount of submerged vegetation can differ. If you had the time/ability to collect wetted width measurements at each interval this could be used to control for this variation by adding it as a covariate in the inferential test. Otherwise, interpret with caution.

Aside from the above, in some cases it can also be difficult to assign proportions consistently in the field, such as when invasive reed canarygras forms dense mats that obscure banks and underlying soils/substrates. If there were such difficulties in the field, these should be clearly stated to aide interpretation. For the purposes of interpretation, remember that in this section we are considering the proportions of the length and the area of residual pools that are occupied by aquatic vegetation. As such, the field observations should relate to proportions of wetted width (not bankfull), because at flows approaching zero the vegetation between observed wetted width and bankfull width would not occupy any part of the residual pool.

🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿 🌿🌿🌿🌿🌿

The tables generated and exported by the code below will show the proportional length of residual pools in which aquatic vegetation was present (regardless of quantity/width), in the same way that we used the fines data. The tables also show the plan area of residual pools that featured aquatic vegetation, based on your estimates of the proportion of wetted width occupied by vegetation at each interval, integrated along the length of each pool.

Just like for fines, for aquatic vegetation we default to using the full contingent of pools. Our rationale is that vegetation (particularly certain problematic invasive species) can be involved in the process of infilling pools by trapping sediment, and it may be of restoration significance to track changes in vegetation in minor nascent or infilling pools. However, as before, we have included in the code below the ability to select a subset of pools, for example if you want to focus on factors such as the importance of aquatic vegetative cover that would plausibly be used by rearing fish in summer (i.e., only in pools of a certain depth). Our default is 0 (no subsetting), and if you change it, ensure that when comparing systems or years you are comparing like for like.

# Define custom threshold if required. We recommend 0 (all residual pools, regardless of maximum depth)
custom_depth <- 0  # Replace with your desired threshold in metres, or set to NULL to use the standard deviation of the reach (that year)

Reachwide Submerged Vegetation for Each Site-Year

	site	year	pool_id	pool_length	veg_length	veg_proportional_length	veg_proportional_area
Brousseau Channel_2024_Total	Brousseau Channel	2024	Total	128	0	0.0000000	0.0000000
Brousseau Channel_2025_Total	Brousseau Channel	2025	Total	142	0	0.0000000	0.0000000
Hatchery Channel_2024_Total	Hatchery Channel	2024	Total	133	112	0.8421053	0.8398496
Hatchery Channel_2025_Total	Hatchery Channel	2025	Total	140	117	0.8357143	0.8321429

Reachwide Emergent Vegetation for Each Site-Year

	site	year	pool_id	pool_length	veg_length	veg_proportional_length	veg_proportional_area
Brousseau Channel_2024_Total	Brousseau Channel	2024	Total	128	0	0	0
Brousseau Channel_2025_Total	Brousseau Channel	2025	Total	142	0	0	0
Hatchery Channel_2024_Total	Hatchery Channel	2024	Total	133	0	0	0
Hatchery Channel_2025_Total	Hatchery Channel	2025	Total	140	0	0	0

Reachwide Floating Vegetation for Each Site-Year

	site	year	pool_id	pool_length	veg_length	veg_proportional_length	veg_proportional_area
Brousseau Channel_2024_Total	Brousseau Channel	2024	Total	128	10	0.0781250	0.0781250
Brousseau Channel_2025_Total	Brousseau Channel	2025	Total	142	10	0.0704225	0.0704225
Hatchery Channel_2024_Total	Hatchery Channel	2024	Total	133	41	0.3082707	0.0748120
Hatchery Channel_2025_Total	Hatchery Channel	2025	Total	140	43	0.3071429	0.0737500

Inferential Tests - Aquatic Veg

Just like with the fines data, we are dealing with data that are bounded by 0 and 1 (for both proportions of length and for area, across each type of vegetation). The non-transformed data would be expected to violate the assumptions of ANOVA, so we omit the initial diagnostic tests. We will again apply Beta regressions.

Below we display the Beta regression model summaries for all vegetation categories (where present), both length and area proportions. We omit the diagnostic plots from the display, but when you run the code for your data, you should examine the plots.

## 
## 
## 
## Model Summary for veg_proportional_length in subveg_result dataset:
## 
## Call:
## betareg(formula = response_adjusted ~ site * Time_Category, data = filtered_data)
## 
## Quantile residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3441 -0.3019 -0.3019  0.4249  0.9525 
## 
## Coefficients (mean model with logit link):
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              -1.594e+00  6.142e-01  -2.596  0.00944
## siteHatchery Channel                      2.766e+00  8.791e-01   3.146  0.00166
## Time_CategoryBefore                      -8.473e-17  6.483e-01   0.000  1.00000
## siteHatchery Channel:Time_CategoryBefore -3.269e-01  9.353e-01  -0.350  0.72668
##                                            
## (Intercept)                              **
## siteHatchery Channel                     **
## Time_CategoryBefore                        
## siteHatchery Channel:Time_CategoryBefore   
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)   0.7401     0.1778   4.162 3.16e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 113.4 on 5 Df
## Pseudo R-squared: 0.7103
## Number of iterations: 22 (BFGS) + 2 (Fisher scoring) 
## 
## 
## 
## Model Summary for veg_proportional_area in subveg_result dataset:
## 
## Call:
## betareg(formula = response_adjusted ~ site * Time_Category, data = filtered_data)
## 
## Quantile residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3905 -0.3036 -0.3036  0.3115  1.0964 
## 
## Coefficients (mean model with logit link):
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              -1.661e+00  6.087e-01  -2.729  0.00635
## siteHatchery Channel                      2.498e+00  8.694e-01   2.873  0.00407
## Time_CategoryBefore                       2.271e-16  6.425e-01   0.000  1.00000
## siteHatchery Channel:Time_CategoryBefore -7.093e-02  9.408e-01  -0.075  0.93990
##                                            
## (Intercept)                              **
## siteHatchery Channel                     **
## Time_CategoryBefore                        
## siteHatchery Channel:Time_CategoryBefore   
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)   0.8011     0.1898    4.22 2.44e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 103.4 on 5 Df
## Pseudo R-squared: 0.6998
## Number of iterations: 27 (BFGS) + 2 (Fisher scoring) 
## 
## Warning: optimization failed to converge for veg_proportional_length in emergveg_result dataset
## 
## 
## 
## No model convergence for veg_proportional_length in emergveg_result 
## 
## Warning: optimization failed to converge for veg_proportional_area in emergveg_result dataset
## 
## 
## 
## No model convergence for veg_proportional_area in emergveg_result 
## 
## 
## 
## Model Summary for veg_proportional_length in floatveg_result dataset:
## 
## Call:
## betareg(formula = response_adjusted ~ site * Time_Category, data = filtered_data)
## 
## Quantile residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1733 -0.5225 -0.3722  0.0608  1.9359 
## 
## Coefficients (mean model with logit link):
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              -1.96714    0.59286  -3.318 0.000907
## siteHatchery Channel                      0.64810    0.74930   0.865 0.387073
## Time_CategoryBefore                      -0.16239    0.61162  -0.266 0.790619
## siteHatchery Channel:Time_CategoryBefore  0.06658    0.87227   0.076 0.939159
##                                             
## (Intercept)                              ***
## siteHatchery Channel                        
## Time_CategoryBefore                         
## siteHatchery Channel:Time_CategoryBefore    
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)   1.5015     0.4167   3.603 0.000314 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 91.44 on 5 Df
## Pseudo R-squared: 0.2478
## Number of iterations: 21 (BFGS) + 1 (Fisher scoring) 
## 
## 
## 
## Model Summary for veg_proportional_area in floatveg_result dataset:
## 
## Call:
## betareg(formula = response_adjusted ~ site * Time_Category, data = filtered_data)
## 
## Quantile residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9762 -0.6076 -0.3915 -0.2771  2.7722 
## 
## Coefficients (mean model with logit link):
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -3.0740     0.5594  -5.495  3.9e-08
## siteHatchery Channel                       0.3284     0.6659   0.493    0.622
## Time_CategoryBefore                       -0.2105     0.5587  -0.377    0.706
## siteHatchery Channel:Time_CategoryBefore   0.1576     0.7767   0.203    0.839
##                                             
## (Intercept)                              ***
## siteHatchery Channel                        
## Time_CategoryBefore                         
## siteHatchery Channel:Time_CategoryBefore    
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)   
## (phi)    6.103      2.047   2.981  0.00288 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 109.5 on 5 Df
## Pseudo R-squared: 0.1604
## Number of iterations: 23 (BFGS) + 2 (Fisher scoring)

Here is that basic interpretation assistance from before, but also remember that there is extra nuance now you are dealing with proportions of pool length or proportion of pool plan area.

• “Optimization failed to converge…” If you see this warning, it may be that there were 0 values across an entire Site. This should be fairly easy to interpret without inferential tests.

• Quantile residuals measure the difference between the actual data and predicted values. The median should be close to 0. Also check the residuals vs fitted plot to ensure there is no clear trend or unusual outliers

• Coefficients (mean model): these are the effects of the predictors (site, time) on the response variable (fine sediment proportion in pools). Look for significant p values in the last column, particularly for the Site_____:Time_Category_____ row, which is the interaction term (if significant, this means the proportion of fines in residual pools changed over time differently between sites).

• Phi coefficients: there is a p-value here as well, but this refers to whether the spread of the data around the mean was an important feature. If the phi estimate is very high, this means the data was concentrated around the mean (fines proportion doesn’t fluctuate much between sites/times). Consider phi of 0.1 to be spread out data, 1 to be moderate, and 10 to be concentrated around the mean.

• Pseudo R-squared: this is how much of the variation in fines in pools can be explained by the site and time variables. Often low in ecological contexts.

Next up are the promised plots that show the distribution of fines and aquatics across the reach, which should facilitate interpretation.

Composite Figure

This figure juxtaposes the plots showing the residual surfaces with the presence/absence record of fine sediments, and with the wetted width proportion of the three categories of aquatic vegetation. Where factors like excess sedimentation, eutrophication, invasive aquatics, excessively fluctuating dissolved oxygen, and/or lack of pool habitat is a concern, this figure should be useful to visually compare reaches or track changes over time.

Remember that in a BACI study our primary interest is whether our restoration site changed in a different way than the control/reference sites changed across time - this is not always easy to see in figures, so consider the figure alongside the inferential test results from above. Bear in mind that the figure presents data by year whereas our inferential tests consider before vs after time categories.

If you have several sites/years, you can modify the code below to choose which figures are plotted in RStudio. If left at NULL site and NULL year, the code will display all site-years in one composite figure. You can alter this code to display one year for one site, all years for one site, or all sites for one year.

Regardless of what you view in RStudio, the code will export all individual site-year plots and the full composite figure as html files to your output folder.

# User-defined site and year for individual plot display
selected_site <- NULL # Replace with the desired "site name" (include the quotation marks) or use NULL for all sites
selected_year <- NULL # Replace with the desired year (no brackets) or use NULL for all years

Now that you have visualisations of the residual pool data, and statistical tests to examine the statistical significance of differences in reach pool characteristics among reaches and years, it is time to consider the biological/morphological significance and the implications for your project.

Next Steps

We can now move on to looking at a few more components of habitat that may be useful to monitor for your project. We will export the residual pool data we generated in this code for use later on (when considering habitat unit classifications). This will be exported to your output folder (which by now, might be a bit busy. If moving files around, remember that the html files need to be in the same folder as their similarly-named dependency folders).

References

Clark, C., Roni, P. and Burgess, S., 2019. Response of juvenile salmonids to large wood placement in Columbia River tributaries. Hydrobiologia, 842(1), pp.173-190.

Geissinger, E.A., Khoo, C.L., Richmond, I.C., Faulkner, S.J. and Schneider, D.C., 2022. A case for beta regression in the natural sciences. Ecosphere, 13(2), p.e3940.

Kaufmann, P.R., Levine, P., Peck, D.V., Robison, E.G. and Seeliger, C., 1999. Quantifying physical habitat in wadeable streams (p. 149). USEPA [National Health and Environmental Effects Research Laboratory, Western Ecology Division].

Lisle, T.E., 1986. Effects of woody debris on anadromous salmonid habitat, Prince of Wales Island, southeast Alaska. North American Journal of Fisheries Management, 6(4), pp.538-550.

Lisle, T.E., 1987. Using” residual depths” to monitor pool depths independently of discharge. Res. Note PSW-RN-394. Berkeley, CA: US Department of Agriculture, Forest Service, Pacific Southwest Forest and Range Experiment Station. 4 p, 394.

Mossop, B. and Bradford, M.J., 2006. Using thalweg profiling to assess and monitor juvenile salmon (Oncorhynchus spp.) habitat in small streams. Canadian Journal of Fisheries and Aquatic Sciences, 63(7), pp.1515-1525.

Robison, E.G., 1998. Reach scale sampling metrics and longitudinal pattern adjustment of small streams. PhD Thesis. Oregon State University.

Stack, W.R., 1988. Factors influencing pool morphology in Oregon coastal streams. MSc Thesis. Oregon State University